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I. 


INTRODUCTION 


For  several  years  the  Ballistics  Research  Labora- 
tories, using  computer  codes,  have  been  investigating  the 
phenomena  associated  with  the  penetration  by  a shaped- 
charged  jet  into  vehicle  armor.  The  purpose  of  these  in- 
vestigations has  been  to  provide  an  armor  configuration  that 
is  light  in  weight  but  effective  against  shaped-charged 
jets . 

The  computer  codes  presently  available  can  numeri- 
cally solve  problems  of  normal  or  zero  obliquity  impact. 

For  plane  strain,  planes  of  non-zero  obliquity  can  be 
solved  using  X-Y  versions  of  Eulerian  codes.  There  are 
certain  unique  phenomena  which  occur  during  oblique  impact;  conse- 
quently, a three-dimensional  model  of  the  penetration  pro- 
cess is  required  to  investigate  in  detail  the  phenomena 
associated  with  oblique  impact. 

TRIDORF^  is  capable  of  performing  these  oblique  im- 
pact calculations.  TRIDORF  is  a two-material  version  of 
TRIOIL*  with  a r igid-perf ectly  plastic  strength  formulation. 
The  description  of  the  TRIDORF  code  is  reported  in  Section 
II.  A normal  impact  has  been  checked  against  the  DORF ^ 
code  (two-dimensional)  and  is  reported  in  Section  III. 

A series  of  four  oblique  impact  calculations  (copper 
jet  impacting  a steel  plate)  have  been  solved  using  the 
TRIDORF  code  and  are  reported  in  Section  IV.  In  addition, 
three  oblique  impact  calculations  (plexiglass  spheres  im- 
pacting a plexiglass  plate)  performed  for  the  Naval  Research 
Laboratory  are  reported  in  Section  IV. 

A fuel  target  configuration  has  been  calculated  with 
a two-dimensional  model  (DOKr)  transferring  over  to  the 
three-dimensional  model  (TRIDORF)  when  the  stress  wave 
(generated  by  the  rigid  jet)  reaches  the  walls  of  the 
cylinder . 

* Johnson,  W.  E.,  "TRIOIL  - A Three-Dimensional  Version  of 
the  OIL  Code,"  GAMD-7310,  General  Atomic,  June  1967. 

^ Johnson,  W.  E. , 'TRIDORF  - A Two-Material  Version  of  the 
TRIOIL  Code  with  Strength,"  Computer  Code  Consultants, 
to  be  published. 

^ Johnson  W.  E. , "Code  Correlation  Study ,"  Technical  Report  No. 
AFWL-TR-70-144 , April  1971;  and  Johnson,  W.  E.,  "Develop- 
ment and  Application  of  Computer  Programs  to  Hypervelocity 
Impact,"  3SR-749,  Systems,  Science  and  Software,  July  1971. 
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The  results  of  these  calculations  are  described  in  Section  V. 
The  code  modifications,  for  both  DORF  and  TRIDORF,  required 
for  these  tasks  are  described  in  Section  VII. 

A multi-material  normal  impact  calculation  has  been 
solved  using  the  DORF^  code  and  is  reported  in  Section  VI. 

A splitting  (a  method  of  fractional  steps,  where  one  can 
replace  a complicated  multi-dimensional  problem  by  a suc- 
cession of  simpler  one-dimensional  ones)version  of  the  OIL** 
code  and  techniques  to  eliminate  the  conventional  "if"  tests 
from  Eulerian  formulations  are  described  in  Section  VII. 


II-,  DESCRIPTION  OF  NUMERICAL  TECHNIQUES 

The  TRIOIL^  code  has  been  modified  to  treat  two  materi- 
als, and  a rigid-perf ectly  plastic  strength  model  has  been 
included.  A routine  to  move  massless  tracer  particles  has 
also  been  incorporated  into  the  code.  This  version  is  called 
TRIDORF^  and  is  in  the  final  state  of  documentation.  The  basic 
difference  equations  for  the  Lagrangian  phase,  the  convective 
phase  and  the  strength  phase  are  listed  below.  Familiarity 
with  the  OIL'^  family  of  codes  is  assumed. 


■•Johnson,  W.  E.  , "OIL  - A Continuous  Two-Dimensional 
Eulerian  Hydrodynamic  Code,"  GAMD-5580,  General  Atomic, 
June  1965. 
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A.  Difference  Equations  For  Phase  1 (Lagrangian  Phase) 


The  Eulerian  equations  we  wish  to  solve  are  the 
following : 


+ V • (pu)  = 0 (1) 

+ u • Vpu  = -VP  (2) 

^ + V • (pEil)  - -V  • (Pu)  (3) 


Equation  1 is  che  conservation  of  mass  equation,  Equation 
2 is  the  conservation  of  momentum  equation,  and  Equation  3 
is  the  conservation  of  energy  equation. 

The  second  terms  on  the  left  side  of  Equations  1,  2 
and  3 are  temporarily  dropped.  Their  contributions  are 
later  approximated  when  we  move  mass  across  cell  bound- 
aries. Rewriting  Equations  2 and  3 in  cartesian  coordi- 
nates results  in  Equations  4,  5,  6 and  7. 


9u 

IZ 

(4) 

9t 

9x 

9v 

(5) 

9t  “ " 

9y 

9w  _ 

(6) 

9t  “ “ 

9z 

9E 

9t 

9 (Pu) 
9x 

9 (Pv) 
9y 

9 (Pw) 
9z 

(7) 
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p = density  of  cell  L,  g/cm^ 

X = X coordinate,  cm 

y = y coordinate,  cm 

z = z coordinate,  cm 

u = X component  of  velocity,  cm/shake 

(1  shake  = 10”®  sec) 

V = y component  of  velocity,  cm/shake 

w = z component  of  velocity,  cm/shake 

P = f(p,Q)  from  the  equation  of  state 
= material  pressure,  jerks/cm® 

(1  jerk  = 10 ‘ ® ergs) 

E = total  specific  energy,  jerks/g 

t = time,  shakes 

Q ■ specific  internal  energy,  jerks/g. 

Rewriting  Equation  7, 

- 9 (Pu)  3 (Pv)  _ 3 (Pw) 

^ 9t  “ 3x  " 9y  " 3z 

or,  substituting  for  the  total  specific  energy, 

3 (Pu)  _ 3 (Pv)  _ 3 (Pw) 
3x  3y  3z 

and  differentiating  , 


3 

3t 


Q + 


1 / ' 
2(" 


+ V + w 


The  p,  u,  V,  w,  P,  E and  Q are  all  cell-centered 
quantities  with  the  cell  index  L,  which  is  defined  as 


L = (J-1)  iMAX  + I + (K-1)  iXMAX 


where  iXMAX  = (iMAX)(jMAX) . 
Also, 

1 

x(I)  = 2^  Ax(i) 

i=l 


J 

y(J)  = 53  Ay(j) 

j = l 

K 

z(K)  = 23  AZ(k) 

k=l 


The  L storage  for  the  cell  quantities  is  indicated 
in  the  figure  below  for  a grid  of  4 4 ^ 4 . 
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The  (~)  refers  to  the  velocities  at  the  completion  of 
Phase  1 (not  av  cycle  n+1,  since  we  have  neglected  the 
transport  terms  at  this  stage) . Here  the  acceleration  of 
cell  L is  seen  to  depend  only  on  the  pressures  in  the 
neighboring  cells  and  not  that  of  L.  Defining  pressures 
at  interfaces, 


PBINd”  = 

P2r"  = L+iXMAX  L 


Substitution  of  these  interface  pressures  into  the  three- 
momentum  equations  results  in 


n 

u^  - u^  = 


At 

_n 


PL^  - PHr"^ 


Ax(I) 


= 


AA 

n 


PBLO’^  - PABOVE*^ 


Ay(J) 


w,  - w,  = 


At 

_n 


PBINd’^  - PZR^ 
Az  (K) 


where  the  usual  tilda  (~)  designates  the  new  velocities 
(not  at  cycle  n+1) , since  we  have  temporarily  dropped  the 
transport  terms. 

The  specific  internal  energy  equation  becomes 


At  L 


n+1/2 

L-iMAX 


V 


n+1/2 

L+iMAX 


n+1/2 

^L-1 


2Ay (J) 


- u 


n+1/2 

L+1 


2Ax(I) 


w 


n+1/2 

L-iXMAX 


- w 


n+1/2 

L+iXMAX 


2Az (K) 


j y^iASS 


MASS  = AMX(L)  + AMD(L) 
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where  Q = specific  internal  energy,  AMX  = X mass  of  cell, 
gm  and  AMD  = • mass  of  cell,  gm.  The  reason  for  the 
velocities  at  time  (n+1/2)  can  be  seen  from  energy  con- 
servation considerations.  Defining 


u 


n+1/2  _ 


v*^  + V 
n+1/2  _ '^L 


w 


n+1/2  _ 


and 


VBLO  = 


'^L-iMAX 


VABOVE  = 


'^L  ^L+iMAX 


UL  = 


+ Vl 


UBIND  = 


w + w 

L-iXMAX 


WZR  = 


W + N 
L L+iXMAX 


then 
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URr"^  + URR  . UBIND*'  + UBIND 
2Ax{I)  2Az  (K) 


WZR  + WZR 
2Az  (K) 


The  solution  of  the  momentum  equations  are  very 
straightforward;  however,  the  solution  to  the  energy 
equation  requires  the  old  and  new  velocities.  Two  passes 
are  made  through  the  routine:  The  first  pass  formulates 

the  interface  velocities  for  the  work  term  contri- 
bution to  the  internal  energy  before  integrating  the 
momentum  equations  (since  we  have  allowed  only  one  array 
per  velocity  component) , and  the  second  pass  bypasses  the 
momentum  equations  and  only  compute  the  new  interface 
velocities  for  their  contribution  to  the  work  term. 

The  change  in  specific  internal  energy  for  each 
material  is  assumed  to  be  proportional  to  the  fractional 
volume  of  each 


AQj^  = AQf/M^ 

and 

AQ.  = AQ(l-f)/M. 

where  (x)  and  (•)  refer  to  the  two  different  materials, 

X and  dot  respectively,  and  f is  the  fraction  of  the  total 
cell  volume  occupied  by  (x)  material.  The  factor,  f,  is 
calculated  from  the  equation  of  state,  where  we  iterate  on 
the  densities  until  the  pressures  of  each  material  are  the 
same . 


B.  Difference  Equations  For  Phase  2 (Convective  Phase) 

The  Eulerian  equations  we  wish  to  solve  are  the 
following; 


Ae  = _7  . 

3t  ^ 

(pu) 

(9) 

9pu  -*• 

-► 

9t  = • 

Vpu 

(10) 

9PE 

— = -V  • 

9t  ^ 

(pEu) 

(11) 

Rewriting  Equation  9,  the  mass  transport  equation  in 
finite  difference  form  results  in 


,n+l 


Pt 


"5-1  “j-i 


At 


_ . ^I-l  '^I-l 

Ay (J)  Ax  (I) 


^K-l  '^^K-l  ~ 
AZ  (K) 


The  velocities  v,  u and  w are  average  type  velocities 
as  described  in  Reference  1. 

Mass,  the  three  components  of  momentum  and  the 
energy  across  all  six  sides  of  the  parallelpiped  are  cal- 
culated. By  conserving  the  x,  y and  z momentum  and  the 
total  energy,  the  new  velocities  are  calculated  and  the 
new  specific  internal  energy  is  the  difference  between 
the  total  energy  in  the  cell  (after  transport)  and  the 
new  kinetic  energy. 


For  the  mass  fluxes,  we  must  determine  (for  a 
mixed  cell)  how  much  of  each  material  to  move.  Four 
possible  situations  concerned  with  two  materials  arise 

(1)  Material  moving  from  a non-mixed  cell  to  a 
mixed  cell.  This  presents  no  difficulty  or 
modification  to  the  transport  terms. 

EXAMPLE  (non-mixed  to  mixed) : 


x and  • materials 
Aj  (cell  interface  area) 

X material 

Mass  flow  is  from  cell  L to  cell  LA. 

AM  = p V A At 
ij  J 

AM^  = AM 
AM,  = 0 

where  p = density  of  x material  in  cell  L. 
(2)  From  a mixed  to  a pure  cell: 


LA 


L 


material 


and  • material 


The  flux  is  positive  from  cell  L to  cell  LA. 
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The  volume  of  material  to  move  from_cell  L to 
cell  LA  is  VAt  (DX  (I)  DZ  (K)  ) where  V is  the 
weighted  velocity  to  use  in  the  mass  flux  equa- 
tion. The  volume  of  x material  in  cell  L is 
VOL,  the  volume  of  the  dot  material  is  (l-f)VOL. 

A.  If  VAttoXd)  DZ(K))  > {l-f)VOl,  let  all  the 
dot  mass  move  out,  that  is, set  Am,  = AMD(L). 

Now  the  remainder  of  the  volume  to  move  will 
be  assigned  to  the  x flux  as  follows:  We  can 

calculate  the  x mass  to  move  from  conserva- 
tion of  partial  volumes, 


VAt(DX(I)  DZ(K))  = VOL^  + VOL, 


AM 

— - + (l-f)VOL. 
^x 


But 


AMX  (L) 

‘^x  f VOL  ' 


thus 

- |v  it  - a-f)  DV(J)) 

B.  If  VAt(DXd)  DZ(K))  < (l-f)VOL  then  the 
X mass  to  move  AMj^  is  set  to  0 and  the  dot 
mass  to  move  AM,  is  given  by 


AM 


V At  = pV  AAt  . 

(l-f)DY(J) 


Material  moving  from  a mixed  cell  to  a mixed 
cell  requires  some  modification  in  the  trans- 
port terms  in  order  to  keep  the  material 
interface  defined  in  a single  mixed  cell. 


EXAMPLE  (mixed  to  mixed) : 


r 


X and 


A 


J 


X and 


• materials 


• materials 


Mass  flow  is  from  cell  L to  cell  LA. 
AM  = pv  AjAt,  where  p is  the  total 
density  of  both  materials  in  cell  L, 
then 


AM  = 


M^  + M^^ 
X X 


^ + M^  + M^ 

A X • • 


AM 


AM  = 


M^  + + M^  + 

XX  . • 


AM 


(4)  From  a mixed  to  a pure  cell,  with  a pure  like 
material  below  - use  the  mass  weighting  scheme 
as  calculated  in  situation  Number  3. 


LA 

L 

LB 


• material 


A 


J 


X and  • material 


* material 
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C.  Difference  Equations  For  Phase  3 (Strength  Phase) 


The  Eulerian  equations  for  the  rigid-perf ectly 
plastic  strength  equations  are  again,  dropping  the  trans- 
port terms : 

Conservation  of  Momentum: 


3u . 

P TT  “ '^ij  , j 


Conservation  of  Energy: 


3E 
» St 


(Olj  U.),  j 


where  E = total  energy  per  gram,  internal  plus  kinetic. 
The  deviatoric  stresses  are  defined  as 


I . . = B (e . . - 6 . . -^1  = B c 
i:  \ i:  1' 


'aa  \ _ 
-D  3 I 


i: 


where 


B = 


y /2 


^ab  ^ab 


\ 


(y  is  defined  on  page  22) 


ID 


3u 

1 /du  3v\ 

1 /3u 

3x 

2 \9y  3x/ 

2^32 

1 /3v  3u\ 

3v 

1 /3v 

2\3x  3y/ 

3y 

2 ^32 

1 /dw  3u\ 

1 /3w  3v\ 

3w 

2 \3x  3z) 

2 Uy  32) 

3z 

9w\ 

3x) 
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-*•  Pm 


2 


+ lx 


1 9u  5v  9w|^ 

3 3x  3y  3z  I 


+ 


I /lx 

'2\3x 


+ 


lx\^ 

dyj 


+ 


i/lx 

2\3z 


+ 


lx  4-  1/lx  -I.  lx  \ 

3x  j 2^3z  3yj 


The  forces  these  stresses  produce  are  then  simply 
the  stresses  multiplied  by  the  area  over  which  they  are 
acting . 


These  stresses  are  indicated  in  the  following 

figure . 
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f 

A 


In  addition  to  the  nine  stresses  as  shown,  there  are 
three  on  the  left,  three  on  the  back  and  three  on  the 
bottom  side  of  the  cell,  a total  of  18  stresses. 


There  are  six  forces 
direction  as  shown. 


F2  2B 


acting  on  the  cell  in  the  y 


F22T  = S22T(Dx (I) ) (Dz (K) ) 
F22B  = -S22B(Dx (I) ) (Dz (K) ) 
F12R  = S12R(Dy (J) ) (Dz (K)  ) 
F12L  = -S12L(Dy (J) ) (Dz (K) ) 
F32F  = S32F(Dx (I) ) (Dy (J) ) 
F32B  = -S32B(Dx (I) ) (Dy  (J) ) 


six  forces  act  on  the  cell  in  the  x direction  as 


shown. 


FUR  = SllR(Dy  (J)  ) (Dz  (K)  ) 
FllL  = - SllL(Dy (J) ) (Dz (K) ) 
F21T  = S21T(Dx (I) ) (Dz (K) ) 
F21B  = -S21B(Dx (I) ) (Dz  (K) ) 
F31F  = S31F(Dx(I) )(Dy(J)  ) 
F31B  = -S3lB(Dx  (I)  ) (Dy  (J)  ) 


Six  forces  act  on  the 

shown . 


cell  in  the  z direction  as 

F33F  = S33F(Dx (I) ) (Dy (J) ) 

F33B  = -S33B(Dx(I) ) (Dy (J) ) 

F13R  = S13R(Dy (J) ) (Dz (K) ) 

F13L  = -S13L(Dy (J) ) (Dz fK) ) 
F23T  = S23T(Dx (i:  ) (Dz  (K)  ) 

F23B  = -S23B(DX (I) ) (Dz (K) ) 


The  six  forces  in  the  z direction  produce  an  ac- 
celeration of 


dw 

dt 


n-i—  [F33F  + F33B  + F13R  + R13L  + F23T  + F23B] . 
''(W 


The  six  forces  in  the  x direction  produce  an  ac- 
celeration of 


du  _ _ 

at  M 


1 

(L) 


(FUR  + FllL  + F21T  + F21B  + F31F  + F31B] 
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and  the  six  forces  in  the  y direction  produce  an  accelera- 
tion of 


dv 
3t  " 


[F22T  + F22B  + F12R  + F12L  + F32F  + F32B]  . 


These  18  forces  will  do  work  such  that 


or 


3E  , 

dt  1]  1 


M 

(L)  9t 


'(L)  * 5 ‘“(M  * ''(L)  * "(L) 


= {FllR(UR)  + FllL(ULL)  + F21T(UA)  + F21B(UB) 


+ F31F(UF)  + F31B(UBA)  + F22T(VA) 


+ F22B(VB)  + F12R(VR)  + F12L(VLL) 


+ F32F{VF)  + F32B(VBA)  + F33F(WF) 


+ F33B(WBA)  + F13R(WR)  + F13L(WLL) 


+ F23T(WA)  + F23B(WB)}  = AE 


or 


AQ 


At 


(L)  M 


(L) 


(AE)  - 


* ^^(L)  ^ *^(L)  ^ 
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where  the  u,  v and  w are  Phase  1 velocities  and  the  velocity 
multiplying  the  stress  is  an  average  velocity  (Phase  1) 
across  that  particular  face  at  which  the  stress  is  applied. 
This  requires  the  saving  of  columns  and  slabs  of  u,  v and 
w velocities  from  Phase  1. 

One  requires  only  to  solve  for  the  stresses  at  the 
top,  right  and  front  of  each  cell,  since  the  stresses  have 
been  calculated  at  the  bottom  (from  the  previous  row),  at 
the  left  (from  the  previous  column  sweep)  and  at  the  back 
(from  the  previous  slab  sweep). 

The  stresses  at  the  right  (R)  side  of  the  cell  (see 
the  figure  below)  are 


a . . 


(e.. 


6 . 

1 ■ 


aa 


where 


(y  is  defined  on  page  22) 


On  the  right  side  (R) , the  stresses  are 
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On  the  top  side  (T) , the  stresses  are 


S22T 


1 jdu  . ^ , 9w\[ 
3 \3x  3y  3z/j 


S21T 


§1  lx  + 

2 3x  3y 


S23T 


BT  3v  3w 

2 3z  3y 


On  the  front  side  (F) , the  stresses  are 


S33F 


BF 


3w 

3z 


1 / 3u  3w  \! 

3 \ 3x  3y  3z  ), 


; S32F 

i 


BF  3w  3v 

X 3y  3z 


S31F 


BF  3w  3ju 

2 3x  3z 


The  velocity  qradients  and  the  second  strain  rate 
tensor  Cgb  are  all  calculated  at  ^iie  three  positions 

indicated  and  the  velocities  are  the  tilda  ( ) Phase  1 
velocities . 

Provisions  exist  for  subcycling  this  strength  rou- 
tine. The  At  for  each  subcycle  is  not  Atp^/N  (where  At[j  is 
the  hydrodynamic  time  step)  but  is  calculated  as  follows: 


At 


(N2+1-N) 2 
H N2(N2+1) 
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where  N2  - number  of  subcycles  and  N is  the  particular 
subcycle . 


The  yield  strength  in  simple  shear  is  considered  to 
oe  a tunction  of  pressure  and  energy  and  is  calculated  as 

y (Yq  + ciP)  (1  - Q/Q^) 


where 


Yq  - yield  strength,  jerks/cm^ 
a = constant 

~ energy  to  melt,  jerks/g  . 


The  rigid  plastic  equations 
those  zones  that  have  a p/p  > l.o 
Q ^ Qj  (Qj  is  an  input  number). 


are  only  applied  to 
or  if  p/pp  < 1.0  and 


Those  cells 
changed  by  -^-5%,  are 

(Phase  3) . 


in  which  the  original  velocity  has  not 
bypassed  for  the  strength  calculation 


III.  TWO  THREE-DIMENSIONAL  IMPACT  CALCULATIONS 


A comparison  between  a normal  (two-dimensional  R-Z) 
calculation  and  a normal  (x-y-z)  calculation  has  been  com- 
pleted for  a normal  impact. 

An  infinite  copper  jet  moving  at  a velocity  of 
7.5  X 10^  cm/sec  impacted  an  aluminum  target  of  finite 
thickness.  The  radius  of  the  jet  (for  the  two-dimensional 
calculation)  was  .115  cm,  and  keeping  the  same  cross-sec- 
tional area  for  the  projectile,  resulted  in  a square  cross- 
sectional  rod  of  .204  cm  on  each  side  (Figure  III-l). 

The  Tillotson  form  of  the  equation  of  state  was  used 
both  in  DORF  and  TRIDORF.  The  two-dimensional  calculation 
was  performed  using  a grid  of  45  zones  in  the  radial,  and 
48  zones  in  the  axial  direction.  The  three-dimensional 
calculation  was  performed  using  a grid  of  12  zones  in  the 
X direction,  23  zones  in  the  y direction  and  12  zones  in 
the  z direction.  Note,  it  is  only  necessary  to  solve  one- 
fourth  of  the  problem  in  the  three-dimensional  mode,  because 
of  the  two  planes  of  symmetry. 

The  results  are  shown  (Figure  III-2)  for  the  pressure 
along  the  axis  of  symmetry  for  the  two  calculations  at  a 
time  of  about  .54  usee.  The  agreement  is  very  good,  con- 
sidering the  very  coarse  zoning  in  the  three-dimensional 
grid.  Results  such  as  these  add  to  the  credibility  of  the 
TRIDORF  code. 


IV.  OBLIQUE  IMPACT 

Two  calculations,  using  TRIOIL  (with  strength),  were 
completed  several  years  ago  for  the  Naval  Research  Labora- 
tory. The  configurations  were  as  follows:  A plexiglass  ^ 

sphere  (R  = .3175  cm)  impacting  at  a velocity  of  6.1  x 10^ 
cm/sec  on  a thin  plexiglass  target  of  2.54  cm  thickness. 
The  spherical  projectile  impacted  normal  (Figure  IV-1)  and 
at  45°  to  the  target  (Figure  IV-2)  where  the  angle  is 
measured  between  the  flight  path  and  the  target. 

The  figures  shown  are  for  the  plane  of  symmetry  as 
a function  of  time.  The  sequence  begins  at  the  upper  left 
and  continues  downward.  The  problem  (normal)  was  rezoned 
(as  indicated  on  the  t = 1.68  usee  figure) . 
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A series  of  four  oblique  impacts  have  been  performed 
using  the  TRIDORF  code.  The  four  calculations  involved  a 
copper  jet  (R  = .1  cm)  impacting  a thin  steel  target  of 
1.27  cm  thickness.  The  copper  jet  was  treated  as  a copper 
rod  with  a rectangular  cross  section  (.1772  cm  on  each  side) 
The  impact  velocity  in  all  cases  was  chosen  as  7.55  x 10^ 
cm/ sec . 


Massless  tracer  particles  were  placed  in  an  orthogona 
grid  within  the  steel  plate  and  copper  rod.  The  configura- 
tions at  various  times  and  for  various  planes  are  shown  in 
Figures  IV-3,  IV-4,  IV-5  and  IV-6. 

Figure  IV-3  is  a display  of  the  interaction  for  the 
30®  impact  problem.  The  cross-sectional  display  begins  in 
the  upper  left,  proceeding  downwards,  for  the  Z planes 
(coming  out  of  the  paper) . The  jaggedness  of  the  top  sur- 
face of  the  target  is  due  to  the  large  zones  used  in  the 
calculation.  Fine  zones  were  placed  in  all  three  direc- 
tions near  the  original  target-projectile  interface. 

Figure  IV-4  is  a similar  display  for  the  45°  inter- 
action. Figure  IV-5  displays,  in  addition  to  the  Z values, 
the  interaction  (plane  of  symmetry)  as  a function  of  time. 
The  Z planes  are  displayed  for  every  other  Z value.  These 
results  are  for  the  60°  impact. 

Figure  IV-6  is  a similar  display  for  the  77.5°  im- 
pact, as  a function  of  the  Z planes.  This  calculation  will 
require  considerably  more  computer  time  for  complete  pene- 
tration (if  complete  penetration  will  occur) . 


V.  FUEL  TARGET  INTERACTION 


The  problem  of  a copper  jet  traveling  at  a velocity 
of  5.617  X 10^  cm/sec  penetrating  a cylindrical  tank  of  fuel 
(initial  density  = 1 g/cc)  has  been  completed.  A modified 
continuous  boundary  condition  at  the  bottom  of  the  tank 
allowed  the  copper  jet  to  continue  through  the  tank  at  a 
constant  velocity.  The  copper  jet  was  treated  as  a rigid 
body,  while  the  fuel  was  treated  as  compressible,  using  an 
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equation  of  state  developed  by  John  Dienes^.  The  form  of 
the  equation  of  state  is 

P = A((>  + B(})^  , where  4>  = log(p/pQ) 

The  calculation  (two-dimensional)  was  carried  to  a 
time  of  t = 1.84y*sec  when  the  copper  jet  had  penetrated 
the  top  of  the  cylindrical  tank.  Stress  waves  have  re- 
flected from  the  top  of  the  tank,  but  the  stress  wave  has 
not  reached  the  upper  right  corner  of  the  tank  at  this  time. 

Subroutines  ADJUST  and  CUBIT  were  used  to  rezone  the 
two-dimensional  grid  such  that  the  axial  and  radial  zoning 
will  be  identical  to  the  x and  y zoning  for  the  three- 
dimensional  problem.  Figure  V-1  indicates  the  initial 
configurations  for  the  two-  and  three-dimensional  problems. 
The  radial  pressure  at  the  bottom  of  the  cylindrical  tank 
is  displayed  in  Figure  V-2B.  The  axial  pressures  (for  given 
radial  positions)  are  displayed  in  Figure  V-2A.  These  re- 
sults are  at  the  completion  of  the  two-dimensional  calcu- 
lation. 

Figure  V-3A  shows  density  profiles  for  an  early  time 
(2-D)  interaction.  One  can  see  the  cavitation  appearing 
along  the  rigid  rod. 

Representative  pressures  (at  the  completion  of  the 
three-dimensional  calculation)  are  shown  in  Figures  V-3B 
and  V-3C.  Figure  V-3B  shows  pressure  contours  for  the 
front  face  of  the  three-dimensional  configuration.  Figure 
V-3C  shows  the  pressure  contour  for  the  bottom  face  of  the 
box.  Figure  V-4A  presents  pressure  versus  distance  along 
the  bottom  and  top  edge  of  the  front  face  of  the  cube. 

Figure  V-4B  presents  pressure  versus  time  at  the  top  right 
hand  corner  and  bottom  right  hand  corner  of  the  front  edge 
of  the  box.  Additional  pressure  and  velocity  profiles 
are  necessary  to  complete  the  flow  patterns  in  the  corners 
of  the  box. 


Dienes,  John#  "The  Response  of  Gelatin  to  High  Speed 
Impact,  SSS-IR-74-2120. 


VI.  MULTI -MATERIAL  IMPACT 


The  configuration  for  the  multi-material  impact  is 
shown  in  Figure  VI-1. 

A copper  jet  (R  = 1 mm)  impacted  on  a finite  aluminum 
target  (thickness  of  3.75  cm)  at  a velocity  of  7.5  x 10® 
cm/sec.  There  are  35  zones  in  the  radial  direction  and  90 
zones  in  the  axial  direction. 

Figure  VI-2  shows  density  plots  of  the  interaction 
as  a function  of  time.  The  density  plots  are  produced  by 
assigning  'v-lOO  random  particles  to  the  cell  with  the  largest 
density.  All  other  cells  are  then  plotted  with  a number  of 
particles  proportional  to  the  cell  density  divided  by  the 
maximum  density  multiplied  by  100. 
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VII.  CODE  MODIFICATIONS 

Several  code  modifications  were  required  in  order  to 
complete  the  fuel  target  interaction  problem. 

The  slug  was  treated  as  a rigid  body  and  moved  through 
the  grid  with  a constant  velocity.  This  required  modifica- 
tions in  the  transport  phase  and  the  Lagrangian  phase  of 
the  DORF  code.  A continuous  bottom  boundary  condition  was 
required  in  order  to  simulate  the  infinite  slug.  An  arti- 
ficial viscosity  routine  {with  linear  and  quadratic  terms)  was 
incorporated  into  the  DORF  code.  The  viscosity  is  calcu- 
lated at  the  zone  interfaces  rather  than  cell  centers. 

These  artificial  viscosity  terms  were  then  incorporated  into 
the  momentum  and  energy  equations  for  the  Lagrangian  phase 
of  the  calculation. 

The  two-dimensional  calculation  for  the  slug  penetra- 
tion was  carried  out  until  the  stress  wave  reached  the  wall. 

At  this  time,  the  codes  ADJUST^  and  CUBIT ^ were  used  to 
rezone  the  problem  such  that  the  zoning  in  the  R and  Z 
direction  for  the  two-dimensional  results  would  be  identi- 
cal to  the  zoning  in  the  X and  Y direction  for  a plane  in 
the  three-dimensional  grid.  One  can  then  transfer  the  v 
component  of  velocity,  the  specific  internal  energy  and  the 
density  directly  and  transform  the  u component  of  velocity 
from  the  two-dimensional  grid  to  the  u and  w components  in 
the  three-dimensional  grid. 


Modifications  to  the  DORF  code  for  interior  rigid 
boundaries  have  been  completed.  These  interior  boundaries 
must  coincide  with  cell  boundaries.  A material  flag,  "9", 
designates  that  the  4 sides  of  the  cell  are  to  be  treated 
as  reflective  boundaries.  Options  of  allowing  the  bottom, 
top  and  right  side  of  the  computational  grid  to  be  either 
reflective  or  transmittive  have  been  added.  The  sound 
speed  for  the  Tillotson  formulation  of  the  equation  of 
state  has  been  changed  to  perform  the  evaluation  of  the 
analytical  form  of 


c2  ^ 


3P 

9p 


3P 

8Q 


P P 


Basic  equations  for  a splitting  version  of  the  OIL 
code  have  been  derived.  In  addition,  the  conventional  "if" 
tests  have  been  removed  and  replaced  with  functions  whose 
values  are  0 or  1.0.  The  splitting  technique  arid  the 
removal  of  the  conventional  "if"  tests  are  applicable  to 
the  new  generation  of  parallel  and  vector  computers. 
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Figure  IV-1.  Density  as  a Function  of  Time  for  the  Normal 
Impact  (continued) 
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Figure  IV-2.  Density  as  a Function  of  Time  for  the  45° 
Impact 
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Figure  V-1 . 


Initial  Fuel  Target  Interaction  (2D  - 3D) 
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Figure  V-2. 


Pressures  in  3D  Mesh  at  Time  of  Rezone 
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Figure  VI-1.  Initial  Configuration  for  the  Multi-Material 
Normal  Impact 
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